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ABSTRACT 


A  new  formulation  of  electromagnetic  wave  scattering  by  convex, 
two  dimensional  conducting  bodies  was  introduced  by  Kriegsman  et  al 
in  ref.[l].  This  formulation,  called  the  On  Surface  Radiation 
Condition  (OSRC)  approach  is  based  upon  an  expansion  of  the 
radiation  condition  applied  directly  on  the  surface  of  the 
scatterer.  Past  approaches  involved  applying  a  radiation  condition 
at  some  distance  from  the  scatterer  in  order  to  achieve  a  nearly 
reflection  free  truncation  of  a  finite  difference  time  domain 
lattice.  Application  of  a  suitable  radiation  condition  directly  on 
the  surface  of  the  convex  conducting  scatterer  can  lead  to 
substantial  simplification  of  the  frequency  domain  integral 
equation  for  the  scattered  field,  which  is  reduced  to  just  a  line 
integral.  In  this  thesis,  the  application  of  the  second  order 
radiation  boundary  condition  of  ref.[l]  to  concave  and  re-entrant 
two  dimensional  bodies  is  considered.  The  specific  shapes 
considered  are  slit  circular  cylinders  and  thick  semi-circular 
cylinders.  The  formulation  is  done  for  the  more  challenging  case  of 
Transverse  Electric  (TE)  polarization.  The  bodies  considered  here 
are  two  dimensional  in  nature  and  are  perfectly  conducting.  Results 
are  presented  for  the  surface  current  distribution  and  scattering 
width  of  these  cylinders  and  comparison  is  made  with  known  results. 
It  was  found  that  the  OSRC  approach  may  not  be  applicable  to 
concave  structures  and  also  for  higher  values  of  kb  for  convex 
bodies . 
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I . INTRODUCTION 


The  On-Surface  Radiation  Condition  (OSRC)  approach 
presented  here  is  a  medium  frequency  technique  for  modeling 
electromagnetic  scattering  from  complex  bodies.  It  is  based 
upon  an  expansion  of  the  radiation  condition  applied  directly 
on  the  surface  of  the  scatterer.  Past  techniques  involved 
applying  a  radiation  condition  at  some  distance  from  the 
scatterer  in  order  to  achieve  a  nearly  reflection  free 
truncation  of  a  finite  time  domain  lattice.  Application  of  a 
suitable  radiation  condition  directly  on  the  surface  of  the 
scatterer  can  lead  to  substantial  simplification  of  the 
frequency  domain  integral  equation  for  the  scattered  field, 
which  is  reduced  to  just  a  line  integral.  The  OSRC  technique 
converts  the  usual  surface  integral  equation  for  the 
scattering  problem  into  either  an  integration  of  unknown 
quantities  or  a  simple  ordinary  differential  equation  for 
convex  two-dimensional  objects.  For  the  Transverse  Electric 
(TE)  case,  the  problem  reduces  to  solving  an  ordinary 
differential  equation  around  the  scatterer  contour  and  for  the 
Transverse  Magnetic  (TM)  case,  the  solution  is  known 
explicitly.  It  is  currently  applicable  to  convex  conducting 
cylinders  of  arbitrary  cross  sections  for  both  TE  and  TM 
cases.  The  OSRC  approach  was  motivated,  as  mentioned  in 
ref.[l],  by  numerical  experiments  conducted  over  the  past 
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twenty  years  aimed  at  simulating  scalar  or  vector  wave 
propagation  and  scattering  using  a  finite  difference  time 
domain  (FD-TD)  model  of  the  governing  wave  equation. 

In  the  subsequent  sections  of  this  thesis,  we  will  use  the 
OSRC  theory  to  study  scattering  by  two  dimensional,  convex  and 
concave  shaped  conducting  scatterers  for  the  TE  case.  The 
specific  shapes  considered  are  the  slit  circular  cylinder  and 
the  thick  semi-circular  cylinder.  More  emphasis  is  placed  on 
the  TE  case  using  the  second  order  boundary  condition.  Results 
and  comparisons  are  shown  for  the  TE  case.  Comparisons  of 
current  distributions  for  circular  cylinders  with  nomalized 
circumferences  (kb)  of  5  and  10  with  ref.[l]  are  made.  Also 
comparisons  of  bistatic  cross  sections  for  circular  cylinders 
are  made  for  several  values  of  Kb.  For  the  thick  semi-circular 
cylinder  and  the  slit  cylinders,  comparison  is  shown  for  the 
induced  current  density  and  the  normalized  scattering  width. 
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II.  SCATTERING  FROM  CIRCULAR  CYLINDERS:  OSRC  APPROACH 


In  this  chapter  we  will  first  apply  the  OSRC  technique  to 
two  dimensional,  perfectly  conducting  and  convex  shaped 
cylinders.  Next  we  will  apply  the  second  order  boundary 
conditions  to  circular  cylinders  with  radius  r^^  and  obtain  a 
second  order  differential  equation  for  the  surface  scattered 
field.  The  differential  equation  is  solved  by  the  Fourier 
series  method.  Finally  expressions  are  obtained  for  the 
surface  current  distribution  and  the  scattered  fields. 

A.  GENERAL  FORMULATION  FOR  THE  TE  POLARIZATION  CASE 

This  formulation  is  based  on  the  second  order  radiation 
conditions  mentioned  in  ref.[l].  For  the  TE  polerization  case, 
the  problem  reduces  to  solving  an  ordinary  differential 
equation  around  the  scatterer  surface  contour. 

Consider  a  plane  electromagnetic  wave  illuminating  a  two 
dimensional,  perfectly  conducting,  convex  shaped  cylinder  for 
the  TE  polarization  case.  The  incident  wave,  propagating  at  an 
angle  a  with  respect  to  the  -x  axis  (see  Fig.  1) ,  is  given  by 

ff  =[/  •  [/  (1) 

"inc  *  »  ^/nc  '  ' 
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Fig.l  TE  incident  field 

where  z  is  the  unit  vector  parallel  to  the  cylinder  axis;  o  is 
the  frequency  of  the  incident  wave  (k=  u/c) ;  c  is  the  is  the 
speed  of  light  in  free  space;  x  and  y  are  the  corresponding 
cartesian  coordinates  in  the  plane  orthogonal  to  2. 

Since  the  problem  is  two  dimensional,  the  scattered 


magnetic  field  (see  Fig. 2)  is  assumed  to  be 


(2) 


The  scattered  field  at  any  distance  r  is  known  from  the 
values  of  the  field  U^(T')  and  its  normal  derivative 
dUsCr'J/dv  *  on  the  surface  of  the  scatterer  as  (ref.  [6]) 

^  dv'  8v' 


y 


Fig. 2  Geometry  of  the  problem 

where  C  represents  the  boundary  of  the  cylinder's  cross 
section;  0/dv'  denotes  an  outward  normal  derivative  on  C;  G  is 
the  two  dimensional  free  space  Green's  function, 
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(4) 


Gir/D 


ikR) 


R=\F-¥‘\=^J {x-x')  *+  iy-y')^  <5) 

T  and  r’  are  co-ordinates  of  the  field  and  source  points,  i.e. 
^=(x,y)  and  P  =  (x' ,  Y' )  »  and  s*  is  the  arc  length  parameter. 

For  the  TE  case,  vanishing  of  the  total  tangential 
electric  field  on  the  surface  of  the  cylinder  produces 

3v'  dv' 

The  objective  is  to  determine  Uj  on  the  surface  of  the 
scatterer.  For  later  purposes,  it  is  necessary  to  obtain  an 
expression  for  the  far-zone  fields.  Using  the  asymptotic  form 
of  the  Hankel  function  it  can  be  shown  that  the  far-zone 
scattered  field  is 


with 


Us 


v/lTcW  c  dv' 


(7) 


(8) 


where 

cos5=^'.f  and  cost|r=  (H. . 

The  true  fields  satisfy  Sommerfeld's  radiation  condition 
(ref. [6])  on  an  infinitely  large  outer  cylinder.  However,  in 
the  OSRC  technique,  a  higher  order  boundary  condition  such  as 
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B2  (defined  below)  is  applied  directly  on  the  surface  of  the 
object.  Implicit  in  this  approximation  is  the  fact  that  only 
outgoing  waves  are  assumed  to  exist  on  the  surface.  Thus,  the 
approximation  is  expected  to  have  serious  drawbacks  in 
geometries  where  significant  stored  energy  can  exist  such  as 
in  re-entrant  and  slit  cylinders.  In  the  sections  to  follow, 
we  will  be  dealing  precisely  with  these  situations. 

The  form  of  second  order  boundary  condition  Bj  is  given 
from  (ref.[l])  with  -jk  in  place  of  jk  as 

c  (s^)  ,  C(s')  ,,  1  y_a  .9. 

2  8[jic+C(s')l  2[jk+C(sO]  dv' 

where  C(s')  is  the  curvature  of  the  surface  at  s'. 

The  condition  B2U5  =  0  is  applied  to  circularly  shaped 
cylinders  in  the  next  section. 

B.  APPLICATION  TO  CIRCULAR  CYLINDERS 

At  the  surface  of  the  perfect  electric  conductor,  the 
boundary  conditions  on  the  total  fields  state  that 

iixH=Js  ;  nxE=0  (10) 

where  Jj  is  the  total  surface  current.  For  the  TE  case  the 
scattered  fields  are 
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(11) 


h:=u^ 


and 


1  dUs 

JUHo  dv' 


Hence  we  have  9  (U,^+U5)/3v '=0  and  J^=-(U5+U,^)  on  the 
surface,  where  v'  is  positive  for  convex  surfaces  and  negative 
for  concave  surfaces. 

For  a=0  in  (1),  the  normal  derivative  of  Uj  can  be 
rewritten  as 

^  ~-h  =±jkcos^'U^„^  .  (12) 

9v'  (±)dr'  dr' 


Substituting  the  expression  for  and  (12)  in  (9)  ,  setting 
the  curvature  C(s')  equal  to  1/r^,  and  replacing  by 


gives 


where 


(13) 


C,  =  ±- 


Sk^rl-Bjkr,  and 


Bk^rl-12jkZg-3 


C,=- 


*  8ic*ro-12j)cr„-3 


(14) 


In  C, ,  the  upper  sign  is  valid  for  convex  surfaces  and  the 
lower  one  for  concave  surfaces.  The  above  equation  is  valid  on 
the  surface  r*=r^  with  being  2n  periodic.  Equation  (13)  is 
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linear  with  constant  coefficients  and  can  be  solved  using 
standard  analytical  or  numerical  methods.  Analytical  methods 
such  as  the  method  of  undetermined  coefficients  and  variation 
of  parameters  yielded  complicated  integrals  for  the  solution 
of  the  second  order  non-homogeneous  differential  equation 
(13)  .  Therefore  (13)  was  solved  using  a  Fourier  series 
expansion. 

1.  SOLUTION  OF  THE  DIFFERENTIAL  EQUATION  FOR  THE 
SURFACE  SCATTERED  FIELD 

Rewriting  (11)  in  terms  of  t/g  gives 

^2-^-t^s=C'iCOs4)^t7,„^  • 

The  incident  field  on  the  surface  of  the  cylinder  is 


^inclr'.r, 


_g-iJcr'co8+' 


(16) 


which  may  be  expressed  in  a  Fourier  series  as  (ref. [6]) 

■  ‘I’* 

We  will  assume  a  complex  exponential  form  for  the  Fourier 
series  representation  of  U^, 
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(18) 


Us=  E 


a— « 


The  coefficient  will  be  solved  by  substituting  into 
equation  (15) .  For  Ug  on  the  circular  portion  of  the  cylinder, 
we  have 


dU. 


(19) 


and 


a^-  " 


(20) 


Substituting  these  into  the  Left  Hand  Side  (LHS)  of  (15) , 
we  get 


X2«.« 


(21) 


Now  expanding  the  Right  Hand  Side  (RHS)  of  (15) , 


C,U,„,cos4>'=cXj-JJkr,)  (e^’+'+e'^*') 


We  now  let  n+l=in  -»  n=in-l  in  the  first  series  and  n-l=in 
n=in+l  in  the  second  to  obtain 
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(23) 


In  obtaining  (23)  ,  use  was  made  of  the  expression 
.  Equating  the  coefficients  of  (21)  and  (23) 


yields  the  following  expression  for  the  Fourier  series 
coefficient  A^, 


l  +  Cjil* 


(24) 


We  note  that  the  coefficients  are  dependent  on  the  radius  r^, 
of  the  cylindrical  surface.  For  numerical  calculations,  the 
series  in  (18)  needs  to  be  truncated.  For  determining  the 
truncation  value  of  n,  we  determine  A„  for  large  n.  For  large 
values  of  the  order  n,  we  have  from  ref. flO] 

Jjkr„)  -  _i:_(.f^)".  (25) 

2n 


Also  using  the  recurrence  relation 


-*4.1  1  ' 


(26) 


we  obtain  for 


ekr^ 

~2ir 


<  1 
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(27) 


2^2icn  2/j 


(lh£o)n-l 


Therefore  for  large  n,  can  be  expressed  as 


2,/2x  2 


n 


(28) 


Hence  decays  very  rapidly,  for  example 


!io  _  3^0-11 


2.  SURFACE  CURRENT  DISTRIBUTION 

Having  found  the  surface  magnetic  field,  it  is  a  simple 
task  to  determine  the  surface  current.  The  total  electric 
surface  current  for  the  circular  cylinder  is  given  by 

•  (29) 

Substituting  (17)  and  (18)  in  (29)  gives 

«• 

.  (30) 

After  expanding  and  simplifying,  (30)  reduces  to 
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{krj  1  cosiKb' 


(31) 


where  e„  =  1  for  n  equal  to  zero  and  c„  =  2  for  n  not  equal  to 
zero,  and  C,  and  Cj  are  given  in  (14)  .  However  for  the 
positive  sign  is  taken  since  the  surface  is  convex.  The 
magnitude  of  the  normalized  surface  current  can  be  expressed 
as 


-7^  1=  |-£  i  1  COSn^'\  . 

tilne  0  l  + 


3.  FAR  ZONE  SCATTERED  FIELDS 

The  far-zone  scattered  magnetic  field  for  the  TE 
case  given  at  (8)  can  also  be  expressed  as  (ref. [5]  page  374) 

^  ^  4  nk  Jc  dv' 

where  C  is  a  closed  contour  around  the  cylinder. 

In  our  case,  with  reference  to  Fig. 2,  we  have 
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9 


(34) 


g  Jkt.  r^_  ^Jkz„coa  ( V-*) 


(35) 


(jp.  O')  =±COS  (<j>'-(J))  , 


(36) 


0//f 

av' 


(37) 


Substituting  (34)  through  (37)  in  (33)  and  simplifying  we 

In  the  derivation  of  the  above  equation,  we  made  use  of 

the  fact  that  yJlcos((}.'-<l»)«'^-*“^*''*>»-i-U,(<^-“^*'**’) . 

dr'  * 


For  a  full  cylinder  we  can  further  simplify  equation  (38) 
by  solving  the  integral  portion  as 

2  R  —JnM 

/e''*^-”’(*'-*’e^'^d(t)'=(-l)"2iie  »  .  (39) 

0 

Substituting  the  above  value  in  (38)  and  after  some 
manipulation,  we  get 
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(40) 


+j  ^A„y„{kr^)]co3n^  . 


We  list  below  a  summary  of  the  various  quantities  that 
were  determined  in  this  chapter. 


Surface  scattered  field.  .  .U, 


/=r  -■ 


Slr'=r, 


Total  surface  field.  .  •  1 

•• 

Surface  current . ‘^♦“"52  [Aj^+j'^Jj^ikr^)]  e^"*' 


Fourier  coefficient . A„= 


CJ-(^*^^ji^(kr,) 


s  ,  e'i’^ 


-  ■.•-(«♦!)  1  /  1  \  2 
C,  2v^  2n  '  'n 

llrroCoa  (♦'-♦) 


Far  fields.  . 


/r 


In  the  next  chapter  we  will  develop  expressions  for  the 
case  when  the  cylinder  is  comprised  of  more  than  one  circular 
surface. 
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III.  SCATTERING  FROM  CIRCULAR  SLIT  CYLINDERS:  OSRC  APPROACH 


In  this  chapter  we  will  consider  the  case  of  a  slit 
cylinder  with  finite  wall  thickness.  It  is  assximed  that  the 
thickness  of  the  wall  is  electrically  small  so  that  the  fields 
do  not  vary  along  it.  Furthermore,  we  ignore  the  effect  of  the 
edges.  Fig. 3  shows  the  geometry  of  the  slit  circular  cylinder. 
On  the  outer  surface  C^,  the  body  is  convex,  while  on  the 
inner  surface  Cj  it  is  concave.  Since  both  outward  and  inward 
travelling  waves  can  exist  in  the  interior  region,  and  since 
the  OSRC  technique  assumes  waves  of  only  one  kind,  it  is  not 
expected  to  yield  correct  results  for  the  geometry.  We  will 
use  the  results  developed  in  the  previous  sections  for  the 
circular  surfaces.  To  distinguish  between  the  outer  and  inner 
surface  we  shall  use  superscripts  (b  for  the  outer  radius  and 
a  for  the  inner  radius)  for  the  Fourier  coefficients. 

A.  DETERMINATION  OF  SURFACE  CURRENTS 

The  expression  for  the  current  distribution  is  given  in 
(29)  i.e.  =  -(Uj  +  Ujnc) '  where 
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Fig. 3  Circular  slit  cylinder 


In  Fig.  3,  0p,  and  Bpj  are  the  first  and  second  slit  angles 
from  the  -x  axis  respectively;  b  and  a  are  the  outer  and  inner 
radii  of  the  cylinder  respectively;  s  is  the  arc  length  and  is 
taken  as  zero  at  0'=Jr-0p2  at  the  outer  contour  C, ;  the  angles 
a  and  p  are  considered  zero  at  ^'~»r-0p2  and  0'=w^-0p, 
respectively;  s=ba  and  s=a/3  are  the  arc  lengths  of  the  outer 
and  inner  surfaces  respectively. 


For  the  outer  contour  C,  we  get  from  fig.  3, 


and  .e-J"”-'-’--’ 


Similarly  for  the  inner  contour  we  get 


,  IJ  =«-lkaco«*'|  (46) 

and  ^inc  ®  ®  •  '  ' 

Therefore  the  current  distribution  for  the  outer  and  inner 
contours  can  be  obtained  by  substituting  (43)  and  (44)  for  the 
outer  surface,  and  (45)  and  (46)  for  the  inner  surface 
respectively  in  (29) .  After  simplification  the  following 
expression  is  obtained  for  the  magnitude  of  the  normalized 
surface  currents. 

For  the  outer  contour  C,  we  have 

"inc  0 

where  a  varies  from  0  to  2 »r- 0^,-0^,;  c/*”  and  in  are 

given  by  the  expressions  from  (16)  as 


ek^b^-Bjkb  ^{b) _ 4 _  j.Q. 

‘  ”  Bk^b^-12jkb-3  *  "  Bk^b^-12jkb-2 


18 


For  the  inner  contour  Cj  we  have 

^inc  0 

where  /?  varies  from  0  to  STr-e^^-e^, ;  and  are  similar 

except  that  a  negative  sign  is  taken  for  C,  in  (12)  since  the 
inner  contour  is  a  concave  surface. 

Along  the  linear  portions  Cj  and  C^,  we  assume  that  the 
current  varies  linearly  between  the  inner  and  outer  curved 
surfaces . 

B.  DETERMINATION  OF  THE  FAR  FIELDS 

From  Fig. 4  the  total  scattered  magnetic  field  for  the 
hollow  circular  slit  cylinder  with  finite  thickness  will  be 


H 


+ 


♦ 


(50) 


Where  C,  and  Cj  are  the  outer  and  inner  contours  with  radii  b 
and  a  respectively;  0^,-0, j  is  the  total  slit  extent. 
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Y 


Fig. 4  Slit  cylinder  geometry 


From  Fig. 4  the  integration  limits  in  (50)  would  be 


/c=  4  *  4  *  4  *  4 


(51) 


r^^b-^a  ro=a  r^^a^b 
For  the  convex  contour  C,,  (38)  is  used  with  rp=b  and 


A  =A  where 

n  n 


,  (W_ 


cWj-(n*l)j^(jCib) 
l  +  C/*’j3* 


(52) 


Substituting  these  values  in  (38)  we  get 
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(53) 


h!  1  (11)  f  <♦'-♦> 

/r  4  nx  "  Jci 

.ei-fd^/UrA^^U^if  eJtr'co. <*'-*> . 

X  OX'  •'Ci 

The  integral  over  is  evaluated  in  the  following 

fashion: 


f  g  JJdbcos  (♦'-♦) 

Jci 

2n  *“®pl 

_J‘gJlcbcos  (♦'-♦)  gJn*'-  J'  eJ**co8  (♦'-♦)  _ 


(54) 


The  first  integral  can  be  evaluated  in  a  closed  form  as 
in  the  previous  chapter.  In  the  second  integral  we  let 
u=7r-0'-*d0'=-du.  After  some  manipulation  we  get 


r  gjkbcoB  (♦'-♦) 


Substituting  (55)  in  (53)  and  carrying  out  some  manipulation 
we  get 
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(56) 


Hz  =-k 


D  ,2j,\  - 


^  4  ^b['^2ne^iji{kb)Jnikb)*j’^Ai, 


(« 


.Jnikb)'\co3n^-Y^e^**jii{kh)  f  cos  (jiu)  du 

e“ 

.  ®x>i 

+j]^  ”/ (u+<|>)  e"-^*^®"‘“**’cos  (nu)  du]  . 


‘’pi 


For  the  inner  contour  Cj  with  radius  a,  which  has  a 
concave  surface  we  have 


dH. 


dv 


Hl=‘YlAja)e^^  . 
.(a)  Cl^^j-^‘^*^>j'jka) 

(/.  O')  = -cos  ((!>'-({>)  , 
and  f  ds'=f  adi>'=-af  di>'  . 

•  Cj  ^  Cj  V  j 

Substituting  these  values  in  (33)  gives 


//;  =-A:  ® 

‘ci 


jkz „COB  (V-*) 


.  e ^"♦'d4)'+  ^  Ir'-a  < /  e^"*'d(|)')  ] 


(57) 


k^  "ar' 

Following  the  same  procedure  as  for  the  contour  C, ,  we  can 
simplify  the  integral  portion  and  after  some  manipulations  we 
get 
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-jkt  'o'—  ** 

//4=-;c-^  ^  )  2  a  2nc^  [  (-1)  "Jj;  (/ca)  (ka)  -j 

w  ®pi* 

.j'„{ka)  1  co3n(|)-52  ”ji(/ca)  f  e‘-^**®®“*“***cos  (iiu)  du 

e° 

«•  ®Pi 

(-1)”  fcos  (u+(J))  cos  (nu)  du]  . 


(58) 


For  solving  the  linear  portions  and  C^,  we  can  ignore 
the  edge  effects  since  the  OSRC  approach  does  not  currently 
provide  edge  current  singularity  behaviour  as  mentioned  in 
ref.[l].  Therefore,  for  the  linear  portions  we  let  Ug  be  equal 
to  the  average  of  Uj  on  the  contours  C,  and  Cj  at  the  position 
of  the  linear  portion. 

For  contour  Cj  we  have 


Us\c  = 


^5lc/^5lc, 


-J 


^  =  -J-  =-jJtcos  (7i-e_,)  (60) 

dv'  dr' 


and  (jf.  O')  =cos  (<!)'-({))  =cos  (n-0p2-<|))  .  (61) 

Substituting  these  values  in  (33)  and  simplifying  we  get 
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<=-^4(44)  M 


COS0 


P3 


j  (coaepj-coB  (0,^4^^)  ] 


4  nic  cos0pa-cos  (0p,+4>) 

_ g jtMc..e„-c„  (»„.♦,  1  ]  .^  ( A,  *A„  ]  ( .1 , 

—m  2 

[g*J*:«co8(epj4*)_^-iWx;oB(ep,4*)  j  j 


P2 


(62) 


Similarly  for  the  linear  portion  we  have  from  Fig. 4 
U  \  +U  \  "  a  <*>4. a  <«> 


^  =  -jkcos  ( n  -0p, )  e (64) 
(/.O')=cos(n-0p,-<}>)  .  (65) 


Substituting  values  in  (33)  and  after  simplification  we  get 
J.  (^)  ■?  r _ _  r  i>:«lCO80p,-CO8(ep,44)] 

^  4  nJc  ^  cos0pj-cos  (0pj+4>)  ^ 

_gj*McoB0p,-coB(ep,4*)jj  [  A«_::^3e^(_l)«cos720p, 

[g-J*«co8(0p,44)_^-J«>co8(0p,44) j  j  ^  (66) 

In  summary,  the  total  scattered  magnetic  field  for  the 
hollow  circular  slit  cylinder  is 
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_  -  ^itr  -1  o  -f  —  **  , 

"  [-kb['£  2ne„l J„(kb)ji,{kb) 

^X*  ^  0 

+  J  1  cosn4>^-  52e^"ji(Jcb)  r e cos  (nu)  du 

••  ®p* 

(-l)"f  cos  (u+<|>)  cos  (iiu)  6“-^*^®“*“*** du] 

-ica  [  j;  2ne„  [  ( -1)  «cr„  (Ira)  (ka)  -j (ira)  ]  cosn«l> 

0 

e„j  ''Ji,(ka)  f  (nu) du 

°  8pi 

n  ®Pa 

(-!)"[  cos  (u+<t))  e-^**®®*‘“***cos  (/ju)  du] 

"  0pi 

_  r  _ COs9p3 _  IcOS0p,-CO8  <0p,+  0)  1  _  iW>ICOS0p,-CO8  (0pj«^)  J  1 

^  cosOp^-cos  (0p2+4)) 

-  .(*)..  (8) 

-E  e„  ( - 1 ) "  I  "  ,  ]  cosu0p,  (0,..*)  _g-iWx:o. (0p,.*)  3  3 

0  2 

_  r _ COS0pj _  r  iJta(coa0pj-co8(0px+^))_  i*McoBOpi-coB(0p,*4>)l3 

^  COS0p_,-COS  (0pi+(l)) 

(-i)«[.fiL_J!il_]cosn0p,le-^*«°"‘®*>‘*^’-e'^^^^  (67) 

0  2 


The  scattering  width  a^-o  defined  as  (ref.  [6]) 


aa.p=2TCr  ■  — 
l^inc|2 


(68) 


For  the  scattered  field  given  in  (49)  we  have 
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°2-d=-^ \-^  tE  {kb)  Ji,  (kb)  *j  (Ai?)  J  cosn^' 

8pi 

"E  (^)  f  e-^"«»“<“**>cos  (nu)  du 

0  n 

”pt 

••  ®pa 

+jE  (-1)  ”  r  COS  (u+<j))  COS  (JIU)  6“^*^®“  du] 

-ka  2ne„  [  (-1)  ^J„(ka)  J^(ka)  -j-'^Aji^'j^ika)  ]  cosucj) 

0 

"E  ®/v7  f  e-^^*®®“<"*^’cos  (nu)  du 

0  0 

®P1 

+jE  <-1)  ”/ cos  (U+<J))  e^**®®"‘“**>cos  (uu)  du] 

[  J  _  gi*i>[co8ep,-co«  (ep,«^)  j  ^ 


-[ 


cpsB 


“pi 


£2_ 


cos0p2-cos  (0p2+4>) 

-E  Gn  ( - 1 )  ”  t  "  ]  cosiiBp^  t  e - e ]  ] 

n  2 

^£1 _  r  lcosep,-co8  (Opj*^)  1  _  Jict>tcosepi-co8  (e-^**)  1  -I 


COS0 


^  COS0pj-CO3  (0pi  +  <l)) 


•  (69) 
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IV.  NUMERICAL  RESULTS  AND  DISCUSSION 


All  analyses  were  carried  out  for  the  TE  polarization  case 
and  when  the  incident  wave  is  incident  from  the  -x  axis,  i.e. 
from  0=^180”.  In  all  cases  a  and  b  denote  the  inner  and  outer 
radii  of  the  cylinders. 

For  the  circular  cylinder.  Figs.  5  and  6  show  the 
comparisons  of  the  normalized  surface  current  distribution 
between  the  present  formulation  (32)  and  of  ref.[l],  both 
using  the  OSRC  approach  for  kb=5  and  10  respectively.  The 
results  are  comparable,  checking  the  validity  of  our 
formulation.  It  was  shown  in  ref.[l],  that  the  results 
approach  the  exact  solution.  The  computer  program  to  compute 
the  normalized  surface  current  distribution  using  the  OSRC 
approach  is  given  in  Appendix  A.  Figs.  7,8  and  9  show  the 
comparisons  of  the  bistatic  widths  for  the  circular  cylinder 
between  the  present  formulation  and  ref. [8]  for  kb=l,5  and  10 
respectively.  Ref. [8]  uses  the  rigorous  moment  method 
technique  for  determining  the  scattering  cross  sections  of 
infinite  cylinders  of  arbitrary  geometrical  cross  section.  The 
results  again  are  in  good  agreement.  The  computer  program  to 
compute  the  scattering  width  for  circular  cylinders  using  the 
OSRC  approach  is  given  in  Appendix  B.  Figs. 10  and  11  show  the 
comparisons  of  forward  (0*0*)  and  backward  (<^=180  *)  scattering 
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widths  respectively,  versus  kb,  between  the  computed  and 
ref. [5]  (exact  solution).  It  is  seen  that  the  results  due  to 
the  OSRC  approach  show  some  deviation  at  higher  frequencies 
(i.e.  for  large  kb)  in  both  cases.  The  computer  program  to 
compute  a/irh  versus  kb  is  given  in  Appendix  C. 

For  the  hollow  circular  slit  cylinder.  Fig. 12  shows  the 
comparison  of  the  normalized  surface  current  distribution  for 
a  hollow  semi-circular  cylinder  obtained  using  the  present 
formulation  (47)  and  (49)  and  that  in  ref. [7].  The  problem  of 
scattering  in  ref. [7]  is  solved  by  means  of  an  integral 
equation.  In  Fig. 14,  the  current  density,  normalized  to  unity 
for  the  magnetic  field  strength  of  the  incident  wave,  is 
plotted  versus  perimeter  s  normalized  to  unity  for  the  half 
perimeter.  The  distance  variable  s  is  taken  equal  to  zero  at 
the  point  on  the  center  face  nearest  the  source  of  the 
incident  wave.  Moving  clockwise  around  the  perimeter,  s 
increases  in  value  reaching  the  value  s=l  at  the  center  face 
furthest  from  the  source.  The  results  again  are  comparable  to 
the  exact  solution.  The  computer  program  to  compute  (47)  and 
(49)  for  the  circular  slit  cylinder  is  given  at  Appendix  D. 
Fig. 13  shows  the  comparison  of  the  scattering  width.  The 
results  obtained  using  the  OSRC  approach  are  quite  different 
from  the  exact  solutions  as  shown  in  Fig. 13.  This  is  in  spite 
of  the  good  agreement  in  the  magnitude  of  the  current 
distribution.  At  this  point  it  seemed  necessary  to  carry  out 
some  phase  comparisons  for  the  current  distribution,  however 
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no  known  results  were  available.  The  computer  program  in 
Appendix  E  was  thoroughly  checked  for  any  bugs  and  even  tested 
for  circular  cylinders  which  yielded  satisfactory  results.  It 
is  felt  that  the  disagreement  in  the  scattering  width  due  to 
the  OSRC  method  is  problably  due  to  the  incorrect  modeling  of 
phase  of  the  current  distribution.  The  OSRC  technique  assumes 
waves  of  only  one  kind  (outgoing)  to  exist  on  the  scattering 
surfaces  which  may  not  yield  correct  results  for  concave 
structures . 

Fig. 14  shows  comparisons  of  forward  scattering  width 
(0=0")  for  the  hollow  circular  slit  cylinder  (b=a  and  total 
slit  angle  of  10")  between  the  OSRC  approach  and  the  full 
circular  cylinder  from  ref. [5]  (exact  solution).  Fig. 15  shows 
the  comparison  between  the  OSRC  approach  and  ref . [5]  results 
for  the  hollow  circular  slit  cylinder.  Such  shapes  are 
expected  to  have  cross  section  resonances,  indicating  that 
there  is  interior  information  contained  in  the  exterior 
scattering  data.  Here,  again  the  results  are  quite  different. 
Another  peculiar  result  was  that  we  were  getting  the  same 
results  for  forward  and  backward  scattering  widths  using  the 
computer  program  given  in  Appendix  F.  This  might  be  due  to 
some  unforseen  error  in  the  program. 
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- Cai-unded  rcsulls  +♦♦  ref. til  results 


Fig. 5  Comparison  of  surface  current  distribution  on  a  PEC 
circular  cylinder  with  kb=5. 
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Fig.  8  Comparison  of  bistatic  scattering  width  for  the 
circular  cylinder  with  kb=5. 
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Fig.  9  Comparison  of  bistatla  scattering  width  for  the 
circular  cylinder  with  kb=10. 
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- CoMind.ctl  results  ref.tS]  results 
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Fig. 10  Comparison  of  forward  scattering  width  (^=0*)  for  the 
circular  cylinder  between  the  computed  (OSRC  approach)  and 
ref . [5]  (exact  solution). 
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Fig. 11  Comparison  of  backward  scatting  width  (^=180*)  for  the 
circular  cylinder  between  the  computed  (OSRC  approach)  and 
ref. [5]  (exact  solution) . 
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Fig. 13  comparison  of  scattering  width  for  the  hollow  semi¬ 
circular  slit  cylinder. 


Fig. 14  Comparison  of  forward  scattering  width  for  a  hollow 
circular  slit  cylinder  with  bsa  and  total  slit  angle  of  10* 
and  full  circular  cylinder  from  ref . [5]  . 
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-  CoMputed  results  ♦♦+  ref #151  results 


Fig. 15  Comparison  of  forward  scattering  width  for  the  hollow 
circular  slit  cylinder  with  b=a  and  slit  angle  of  10'. 
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V.  CONCLnSIONS 


The  OSRC  approach  was  applied  to  circular  and  hollow  slit 
circular  cylinders.  For  circular  cylinders  the  normalized 
surface  current  distribution  agrees  very  well  with  the  exact 
solution.  Scattering  width  results,  with  different  values  of 
kb  were  also  comparable  with  the  exact  solution,  as  long  as  kb 
was  not  large.  For  the  case  of  hollow  circular  slit  cylinders, 
the  OSRC  normalized  surface  current  distribution  was 
comparable  to  exact  solution  results.  However,  in  case  of  the 
scattering  width,  the  slit  cylinder  results  were  quite 
different  from  the  exact  solution.  This  indicated  a  need  to 
carry  out  phase  comparisons  of  the  OSRC  surface  current 
distribution  with  the  exact  solution. 

It  is  possible  that  the  OSRC  approach  may  not  be 
applicable  to  concave  structures  and  also  for  higher  values  of 
kb  for  convex  bodies.  As  mentioned  earlier,  both  outward  and 
inward  travelling  waves  can  exist  in  the  interior  region  of 
the  concave  scatterer.  Since  the  OSRC  approach  assumes  waves 
of  only  one  kind,  it  is  not  expected  to  be  applicable  to 
concave  bodies.  Further  the  method  fails  when  the  scattering 
surface  has  cavity  shaped  objects. 


APPENDIX  A.  COMPUTER  PROGRAM  TO  COMPUTE  CURRENT  DISTRIBUTION  FOR 
CIRCUZAR  CYLINDERS 


C  THIS  PROGRAM  COMPUTES  THE  SURFACE  CURRENT  DISTRIBUTION  FOR  A 
C  CIRCULAR  CYLINDER  FOR  TRANSVERSE  ELECTRIC  CASE  USING  OSRC 
C  APPROACH. 

REAL*8  J(0:200) ,¥(0:200) ,Dy(0:200) ,DJ(0:200) ,K1D,K2D 
REAL*8  Jl(0:200) ,¥1(0:200) ,DY1(0:200) ,DJ1(0:200) 

REAL  EN,NMX,P,PHI,K,K1,DTR,PI,FR,B,JPHI 
INTEGER  NMAX,N 

COMPLEX  CBl,CB2,i,Sl,S2,S3,S4,DN(999) ,JP 
OPEN ( 3 , FILE= ' DAT . 1 • ) 

WRITE (*,*)  'ENTER  THE  WAVE  NUMBER  K  (REAL)  IN  I/METERS' 
READ(*,*)  K 

WRITE (*,*)  'ENTER  THE  RADIUS  (REAL)  IN  METERS' 

READ(*,*)  B 

WRITE (*,*)  'ENTER  NMAX  (INT) ' 

READ(*,*)  NMAX 
PI=3. 1415927 
DTR=PI/180. 
i=(0,l) 

K1=K*B 

K1D=DBLE(K1) 

CB1= (8 . *Kl**2-8 . *i*Kl)/ (8 . *K1**2-12 . *i*Kl-3 . ) 

CB2=-4 ./ (8 . *K1**2-12 . *i*Kl~3 . ) 

CALL  BES (NMAX, KID, J,Y,DJ,DY) 

C  PLOTTING  ROUTINE 

WRITE (3, 100)  'SURFACE  CURRENT  DISTRIBUTION  FOR  CIRCULAR  C¥L.  ' 

NPHI=91 

PHI 1=0. 

PHI2=180. 

WRITE (3, 110)  NPHI 
WRITE (3, 120)  PHIl 
WRITE (3, 120)  PHI2 
DO  33  P=0. ,180. ,2. 

PHI=DTR*P 

S4=(0,0) 

DO  22  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O. )  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S1=(0,0) 

S2=(0,0) 

S3=(0,0) 

Sl=CBl*i**(-l)  *  DJ(NMX)/(1.+CB2*NMX**2) 

S2=J(NMX) 

S3=COS(NMX*(PI+PHI) ) 

DN(N)=-(EN*i**(-NMX) *(S1+S2) *S3) 
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S4=S4+DN(N) 

22  CONTINUE 

JP=S4 

JPHI=CABS(JP) 
WRITE (3, 120)  JPHI 


33 

CONTINUE 

CLOSE ( 3 ) 

100 

FORMAT (A) 

110 

FORMAT (15) 

120 

FORMAT (E12. 3) 
STOP 

END 

o  o  o  o  o 


APPENDIX  B.  COMPUTER  PROGRAM  TO  COMPUTE  a/b  FOR  CIRCUZAR  CYLINDERS 


THIS  PROGRAM  COMPUTES  THE  SCATTERING  WIDTH  FOR  A  CIRCULAR 
CYLINDER  FOR  THE  TRANSVERSE  ELECTRIC  POLARIZATION  CASE 
USING  THE  ON  -  SURFACE  RADIATION  CONDITION  APPROACH. 
ANGLE  SI  =  PI  -  PHI  GOING  CLOCKWISE  FROM  -180  DEGREES 
ALONG  THE  X-AXIS. 

REAL*8  J(0:200) ,Y(0:200) ,DY(0:200) ,DJ(0:200) ,K1D 
INTEGER  N,NMAX 

REAL  TA,TB,K,K1,K2,TP1,TP2,FIL0N,FR,B,A,L,PI,P,NMX,EN,PHI 
COMPLEX  DN(999) ,S2 ,S3,S4,CN1, i,CBl,CB2 
OPEN(3 , FILE* ' DAT. 2 ' ) 

WRITE (*,*)  'ENTER  THE  WAVE  NUMBER (REAL)  IN  I/METERS' 
READ(*,*)  K 

WRITE (*,*)  'ENTER  THE  OUTER  RADIUS (REAL)  IN  METERS' 
READ(*,*)  B 

WRITE ( * , * )  ' ENTER  NMAX ( INTEGER) ' 

READ(*,*)  NMAX 
PI=3. 1415927 
DTR=PI/180. 
i=(0,l) 

K1=K*B 
K1D*DBLE(K1) 

CB1= (8 . *Kl**2-8 . *i*Kl) / (8 . *K1**2-12 . *i*Kl-3 . ) 

CB2=-4 ./ (8 . *K1**2-12 . *i*Kl-3 . ) 

CALL  BES(NMAX,K1D,J,Y,DJ,DY) 

C  PLOTTING  ROUTINE 

WRITE (3, 960) 'TE  CASE  PEC  CIR.  CYLINDER  SCATTERING  WIDTH  ' 

NPHI=91 

PHI 1=0. 

PHI2=180. 

WRITE (3, 970)  NPHI 
WRITE (3, 980)  PHIl 
WRITE (3, 980)  PHI2 
DO  50  P  =0. ,180. ,2. 

PHI=DTR*P 

C  EVALUATING  THE  FUNCTION  CNl 
S2=(0,0) 

DO  800  N=1,NMAX+1 
NMX=FLOAT (N) -1 . 

S3=(0,0) 

S4=(0,0) 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S3=J(NMX)*DJ(NMX) 

S4=i**(-1)*CB1*DJ(NMX)**2/(1+CB2*NMX**2) 

DN(N) =2 . *PI*EN* (S3+S4) *COS (NMX* (PI-PHI) ) 
S2=S2+DN(N) 
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800  CONTINUE 

CN1=S2 

C  EVALUATING  THE  SCATTERING  WIDTH  DIVIDED  BY  CYLINDER  RADIUS  B 
SIGMA=1./ (4 . *K1) *CABS (-K1*CN1) **2 
WRITE (3, 980)  SIGMA 


50 

CONTINUE 

CLOSE (3) 

960 

FORMAT (A) 

970 

FORMAT (15) 

980 

FORMAT (E12. 3) 

STOP 

END 
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o  o  n  o 


APPENDIX  C.  COMPUTER  PROGRAM  TO  COMPUTE  a/wb  versus  Icb  FOR 
CIRCULAR  CYLINDERS 

THIS  PROGRAM  COMPUTES  THE  SCATTERING  WIDTH  FOR  A  CIRCULAR 
CYLINDER  FOR  THE  TRANSVERSE  ELECTRIC  POLARIZATION  CASE 
USING  THE  ON  -  SURFACE  RADIATION  CONDITION  APPROACH. 
(SIGMA  /PI*B  Vs  KB) 

REAL*8  J(0:200) ,Y(0:200) ,DY(0:200) ,DJ(0;200) ,K1D 
INTEGER  N,NMAX 

REAL  TA,TB,K,K1,K2,TP1,TP2,FIL0N,FR,B,A,L,PI,P,NMX,EN,PHI,R 
COMPLEX  DN(999) , S2 , S3 , S4 , CNl , i, CBl , CB2 
OPEN ( 3 , FI LE= • DAT . 3 ' ) 

WRITE (*,*)  'ENTER  THE  PHI  ANGLE 'REAL)  IN  DECREES' 

READ(*,*)  P 

WRITE (*,*)  ‘ENTER  THE  OUTER  RADIUS (REAL)  IN  METERS' 
READ(*,*)  B 

WRITE (*,*)  'ENTER  NMAX( INTEGER) ' 

READ(*,*)  NMAX 
PI=3. 1415927 
DTR=PI/180. 

PHI  =  DTR*P 
i=(0,l) 

C  PLOTTING  ROUTINE 

WRITE (3, 960)  'TE  CASE  PEC  CIR.  CYLINDER  SCATTERING  WIDTH  ' 

NPHI=59 

PHI 1=0.1 

PHI 2=6. 

WRITE (3, 970)  NPHI 
WRITE (3, 980)  PHIl 
WRITE (3, 980)  PHI2 
DO  50  R  =0.1,6. , .1 
K1=R*B 
K1D=DBLE(K1) 

CB1= (8 . *Kl**2-8 . *i*Kl)/ (8 . *K1**2-12 . *i*Kl-3 . ) 

CB2=-4 . / (8 . *K1**2-12 . *i*Kl-3 . ) 

CALL  BES (NMAX, KID, J,Y,DJ,DY) 

C  EVALUATING  THE  FUNCTION  CNl 
S2=(0,0) 

DO  800  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

IF(NMX.EQ.0.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S3=J (NMX) *DJ (NMX) 

S4=i** (-1 . ) *CB1*DJ (NMX) **2/ ( 1 . +CB2*NMX**2) 

DN(N) =2 . *PI*EN* (S3+S4) *COS (NMX*PHI) 

S2=S2+DN(N) 

800  CONTINUE 
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CN1=S2 

C  EVALUATING  THE  SCATTERING  WIDTH  DIVIDED  BY  CYLINDER  RADIUS  B 
SIGMA«1 ./ (4 . *K1*PI) *CABS (“K1*CN1) **2 
WRITE (3, 980)  SIGMA 
50  CONTINUE 

CLOSE (3) 

960  FOK,'JVT(A) 

970  FORMAT (15) 

980  FORMAT (E12. 3) 

STOP 

END 
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APPENDIX  D.  COMPUTER  PROGRAM  TO  COMPUTE  J,  FOR  SLIT  CIRCULAR 
CYLINDERS 

C  THIS  PROGRAM  COMPUTES  THE  SURFACE  CURRENT  DISTRIBUTION  FOR  A 
C  CIRCULAR  SLIT  CYLINDER  FOR  TRANSVERSE  ELECTRIC  CASE  USING  OSRC 
C  APPROACH. 

REAL*8  J(0:200) ,Y(0:200) ,DJ(0:200) ,DY(0:200) ,K1D,K2D 
REAL*8  Jl(0;200) ,Y1(0:200) ,DJ1(0:200) ,DY1(0:200) 

REAL  EN,NMX,P,PHI,K,K1,K2,SMAX,PI,FR,A,B,JPHI,L,S,S2,SM2,R 
REAL  SMAX2 
INTEGER  NMAX,N 

COMPLEX  CBl,CB2,CAl,CA2,i,Sl,S3,S4,DN(999) ,JP,S5 
OPEN(3 , FILE= ' DAT. 5 ' ) 

WRITE (*,*)  'ENTER  THE  FREQUENCY  (REAL)  IN  MHZ' 

READ(*,*)  FR 

WRITE (*,*)  'ENTER  THE  FIRST  SLIT  ANGLE  (REAL)  IN  DEGREES' 
READ(*,*)  TA 

WRITE (*,*)  'ENTER  THE  SECOND  SLIT  ANGLE  (REAL)  IN  DEGREES' 
READ(*,*)  TB 

WRITE (*,*)  'ENTER  NMAX  (INT) ' 

READ(*,*)  NMAX 
PI=3. 1415927 
i=(0,l) 

DTR=PI/180. 

K=PI*FR/150. 

L=300./FR 
B* . 6^L 
A=.5*L 
K1=K*B 
K2=K*A 
K1D=DBLE(K1) 

K2D=DBLE(K2) 

CALL  BES (NMAX, KID, J,Y,DJ,DY) 

CALL  BES (NMAX, K2D,J1,Y1, DJI, DYl) 

TPl=TA/57.29578 
TP2=TB/57. 29578 
CB1= (8 . *Kl**2-8 . *i*Kl)/ (8 . *K1**2-12 . *i*Kl-3 . ) 

CB2=-4 ./ (8 . *K1**2-12 . *i*Kl-3 . ) 

CA1=- (8 . *K2**2-8 . *i*K2)/ (8 . *K2**2-12 . *i*K2-3 . ) 
CA2=-4./(8.*K2**2-12.*i*K2-3.) 

C  PLOTTING  ROUTINE 

WRITE(3,100)  '  CURRENT  DISTRIBUTION  FOR  CIRCULAR  SLIT  CYL. ' 

NPHI=95 

PHI 1=0. 

PHI2=1. 

WRITE (3, 110)  NPHI 
WRITE(3,120)  PHIl 
WRITE(3,120)  PHI2 

SMAX=B* (2*PI-ABS (TP2-TP1) -PI/2 )+(B-A)+A*(2*PI-ABS(TP2-TPl) -PI/2) 
DO  33  P=0. , .52, .01 
S=P*SMAX/B+PI/2 
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S3=EXP(i*Kl*COS(TP2+S) ) 

S4=(0,0) 

DO  22  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S1=(0,0) 

S2=0. 

Sl=EN*CBl*i**(NMX-l.)  *  DJ(NMX)/(1.+CB2*NMX**2) 
S2=COS(NMX*(S+TP2) ) 

DN(N)=S1*S2 

S4=S4+DN(N) 

22  CONTINUE 

JP=-(S4+S3) 

JPHI=CABS ( JP) 

WRITE (3, 120)  JPHI 
33  CONTINUE 

C  PLOTTING  ROUTINE 

DO  10  R=.57,l. , .01 
S=  (R-.57)*SMAX/A 
S5=EXP(i*K2*COS(TPl-S) ) 

S4=(0,0) 

DO  5  N=1,NMAX+1 

NMX«FLOAT(N)-l. 

IF(NMX.EQ.O. )  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

S1=(0,0) 

S2=0. 

Sl=EN*CAl*i**(NMX-l. ) *DJ1(NMX)/(1.+CA2*NMX**2) 
S2=COS (NMX* (TPl-S) ) 

DN(N)=S1*S2 

S4=S4+DN(N) 

5  CONTINUE 

JP=-(S4+S5) 

JPHI=CABS ( JP) 

WRITE (3, 120)  JPHI 
10  CONTINUE 

CLOSE (3) 

100  FORMAT (A) 

110  FORMAT(I5) 

120  FORMAT (E12. 3) 

STOP 

END 
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APPENDIX  E.  COMPUTER  PROGRAM  FOR  /o/X  versus  0  FOR  SLIT  CIRCULAR 
CYLINDERS 


C  THIS  PROGRAM  COMPUTES  THE  SCATTERING  WIDTH  SQRT( SIGMA/ LAMBDA) 
C  FOR  A  CIRCULAR  HOLLOW  CYLINDER  WITH  A  SLIT  OF  ANY  SIZE  FOR  THE 
C  TRANSVERSE  ELECTRIC  POLARIZATION  CASE  USING  THE  ON  -  SURFACE 
C  RADIATION  CONDITION  APPROACH. 

REAL*8  J(0:200) ,Y(0;200) ,DY(0:200) ,DJ(0:200) ,K1D,K2D 
REAL*8  Jl(0:200) ,Y1(0:200) ,DJ1(0:200) ,DY1(0:200) ,H1 
REAL  EN,NMX,R,U,UU,H,P,PI,PHI,FAC,F(99) ,P1 
INTEGER  M,N,NT,NN,NMAX 

REAL  TA,TB,K,K1,K2,TP1,TP2,FIL0N,FR,B,A,L 

COMPLEX  P2,P3, i,Sl,Tl,T2,T3,T4,CBl,CB2,CAl,CA2 

COMPLEX  C1,C2,C3,C4,C5,C6,C7,C8,CN2,CN3,CN5,CN6,CN7,CN8 

COMPLEX  DN(999) , S2 , S3 , S4 ,S5, CNl , CN4 

EXTERNAL  FI LON 

OPEN ( 16 , FILE= ' DAT . 6 ' ) 

WRITE (*,*)  'ENTER  THE  FREQUENCY (REAL)  IN  MHZ' 

READ(*,*)  FR 

WRITE (*,*)  '  ENTER  THE  FIRST  SLIT  ANGLE (REAL)  IN  DEGREES' 

READ(*,*)  TA 

WRITE (*,*)  'ENTER  SECOND  SLIT  ANGLE (REAL)  IN  DEGREES' 
READ(*,*)  TB 

WRITE (*,*)  'ENTER  NMAX( INTEGER) ' 

READ(*,*)NMAX 
TPl=TA/57. 29578 
TP2=TB/57. 29578 
PI=3. 1415927 
DTR=PI/180. 
i=(0,l) 

K=PI*FR/150. 

L=300./FR 

B=.6*L 

A=.5*L 

K1=K*B 

K1D=DBLE(K1) 

K2=K*A 

K2D=DBLE(K2) 

M=98 

H=ABS (TP2-TP1 ) /FLOAT (M) 

CB1=  (8  .  *Kl**2-8  .  *i*Kl)  /  (8  .  *K1**2-12  .  *i*Kl--3  . ) 
CAl=-(8.*K2**2-8.*i*K2)/(8.*K2**2-12.*i*K2-3.) 

CB2=-4 ./ (8 . *K1**2-12 . *i*Kl-3 . ) 
CA2=-4./(8.*K2**2-12.*i*K2-3.) 

CALL  BES(NMAX,K1D, J,Y,DJ,DY) 

CALL  BES(NMAX,K2D,J1,Y1,DJ1,DY1) 

C  PLOTTING  ROUTINE 

WRITE (16, 960)  'TE  CASE  PEC  SLIT  CYLINDER  SCATTERING  WIDTH  ' 
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NPHI=91 
PHI 1=0. 

PHI2=180. 

WRITE (16, 970)  NPHI 
WRITE (16, 980)  PHIl 
WRITE (16, 980)  PHI2 
DO  50  P  =0.,180.,2. 

PHI=DTR*P 
IF(P.EQ.O.)  THEN 
Tl=i*K* (A-B) *COS (TP2) 

T3=i*K* (A-B) *COS (TPl) 

ELSE 

Tl=COS (TP2)/ (COS (TP2) -COS (TP2+PHI) ) * (EXP(i*K2* (COS (TP2) -COS (TP2+ 
+  PHI) ) )-EXP(i*Kl*(COS(TP2)-COS(TP2+PHI) ) ) ) 

T3=COS (TPl)/ (COS (TPl) -COS (TPl+PHI) ) * (EXP(i*K2* (COS (TPl) -COS (TP1+ 
+  PHI) ) )-EXP(i*Kl* (COS (TPl) -COS (TPl+PHI) ) ) ) 

ENDIF 

T2=EXP( (-i)*K2*COS(TP2+PHI) )-EXP( (-i)*Kl*COS(TP2+PHI) ) 

T4=EXP( (-i) *K2*COS (TPl+PHI) )-EXP( (-i) *K1*C0S (TPl+PHI) ) 

C  EVALUATING  1ST  FUNCTION  F(U) 

DO  100  U  =1,M+1 

UU=TP1+(H*(U-1.)) 

F (U) =COS (Kl*COS (UU+PHI) ) 

CONTINUE 
S1=(0,0) 

P1=0. 

DO  200  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 
EN=2. 

ENDIF 
P2=(0,0) 

PI  =  FILON(F, NMX, TPl, TP2,M+1, 99,1) 

WRITE (*,*)  NMX, PI 
P2=EN*i**NMX*DJ (NMX) 

WRITE (*,*)  P2 
P3=P1*P2 
WRITE(*,*)  P3 
SI  =  SI  +  P3 
CONTINUE 
C1=S1 

C  EVALUATING  2ND  FUNCTION  F(U)  (CN2=Cl-jC2)  &  CN4 


100 


200 
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DO  250  U=1,M+1 

UU=TP1+(H*(U-1.) ) 

F(U)=SIN(K1*C0S(UU+PHI) ) 

250  CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO  300  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

P2=(0,0) 

IF(NMX.EQ.O. )  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

Pl=FILON (F, NMX , TPl , TP2 , M+1 ,99,1) 

S3= ( -1 . ) **NMX*J1 (NMX) *DJ1 (NMX) 

S4=(-l. ) **NMX*i**(-l.) *CA1*DJ1 (NMX) **2/(1. +CA2*NMX**2) 

DN (N) =2 . *PI*EN* (S3-S4 ) *COS (NMX*PHI ) 
P2=EN*i**NMX*DJ (NMX) 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

300  CONTINUE 

CN4=S2 

C2=S1 

CN2*Cl«i*C2 

C  EVALUATING  3RD  FUNCTION  F(U)  &  CN7 

DO  350  U=1,M+1 

UU=TP1+(H*(U-1.)) 

F (U) =COS (UU+PHI) *COS (Kl*COS (UU+PHI) ) 

350  CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO  400  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

S5=(0,0) 

P2=(0,0) 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

Pl=FILON ( F , NMX , TPl , TP2 , M+1 , 99 , 1) 

P2=EN*CBl*i** (-NMX-1. ) * (-1. ) **NMX*DJ (NMX)/ (1 .+CB2*NMX**2) 
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S3=C0S (NMX*TP2) *T2 

S4= (CBl*i** ( -NMX-1 . ) *DJ (NMX) ) / ( 1 . +CB2  *NMX**2 ) 

S5= (CAl*i** ( -NMX-1 . ) *DJ1 (NMX) ) / ( 1 . +CA2*NMX**2 ) 

DN (N) =EN/2 .*(“!.) **NMX* (S4+S5) *S3 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

400  CONTINUE 

C3=S1 

CN7=S2 

C  EVALUATING  4RTH  FUNCTION  (CN3=C3-jC4)  &  CN8 

DO  450  U=1,M+1 

UU=TP1+(H*(U-1.) ) 

F (U) =COS (UU+PHI) *SIN (Kl*COS (UU+PHI) ) 

450  CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO  500  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

S5=(0,0) 

P2=(0,0) 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

Pl=FILON (F, NMX, TP1,TP2,M+1, 99,1) 

P2=EN*CBl*i** ( -NMX-1 .)*(-!.) **NMX*DJ (NMX) / ( 1 . +CB2*NMX**2 ) 
P3=P1*P2 

S3=COS (NMX*TP1) *T4 

S4=(CBl*i**(-NMX-l.)*DJ(NMX) )/(l.+CB2*NMX**2) 
S5=(CAl*i** ( -NMX-1. ) *DJ1 (NMX) )/ (l.+CA2*NMX**2) 

DN (N) =EN/2 . * (-1 . ) **NMX* (S4+S5) *S3 

S1=S1+P3 

S2=S2+DN(N) 

500  CONTINUE 

CN8=S2 
C4=S1 

CN3=C3-i*C4 


C  EVALUATING  5TH  FUNCTION 

DO  550  U=1,M+1 

UU=TP1+(H*(U-1. ) ) 
F(U)=COS(K2*COS (UU+PHI) ) 
550  CONTINUE 

S1=(0,0) 

P1=0. 
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DO  600  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

Pl=FILON ( F , NMX , TPl , TP2 , M+1 , 99 , 1 ) 
P2=EN*i**NMX*DJl (NMX) 

P3=P1*P2 
S1=S1+P3 
600  CONTINUE 

C5=S1 

C  EVALUATING  6TH  FUNCTION  (CN5=C5-jC6) 

DO  650  U=1,M+1 

UU=TP1+(H*(U-1.) ) 

F(U)=SIN(K2*COS(UU+PHI) ) 

650  CONTINUE 

S1=(0,0) 

P1=0. 

DO  700  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

Pl=FILON ( F , NMX , TPl , TP2 , M+1 , 99 , 1) 
P2=EN*i**NMX*DJl (NMX) 

P3=P1*P2 
S1=S1+P3 
700  CONTINUE 

C6=S1 

CN5=C5+i*C6 

C  EVALUATING  THE  7TH  FUNCTION  &  CNl 
DO  750  U=1,M+1 

UU=TP1+(H*(U-1.)) 

F (U) “COS (UU+PHI) *COS (K2*COS (UU+PHI) ) 
750  CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO  800  N=1,NMAX+1 
NMX=FLOAT(N) -1. 

S3=(0,0) 

S4=(0,0) 

IF(NMX.EQ.O. )  THEN 
EN=1. 
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ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

P1=FIL0N ( F , NMX , TPl , TP2 , M+1 , 99 , 1 ) 

P2=EN*CAl*i** ( -NMX-1 .)*(“!•) **NMX*DJ1 (NMX) / ( 1 . +CA2*NMX**2 ) 
S3=J(NMX)*DJ(NMX) 

S4=i**(-1)*CB1*DJ(NMX)**2/(1+CB2*NMX**2) 

DN (N) =2 . *PI*EN* (S3+S4) *COS (NMX*PHI) 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

800  CONTINUE 

C7=S1 
CN1=S2 

C  EVALUATING  THE  8TH  FUNCTION  (CN6=C7-jC8) 

DO  850  U=1,M+1 

UU=TP1+(H*(U-1.) ) 

F (U) =COS (UU+PHI) *SIN(K2*COS (UU+PHI) ) 

850  CONTINUE 

S1=(0,0) 

P1=0. 

DO  900  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

Pl=FILON (F, NMX, TP1,TP2,M+1, 99,1) 

P2=EN*CAl*i** ( -NMX-1. ) * (-1. ) **NMX*DJ1 (NMX) / ( 1 . +CA2*NMX**2 ) 
P3=P1*P2 
S1=S1+P3 
900  CONTINUE 

C8=S1 

CN6=C7+i*C8 

SIGMA=1 ./ (8 . *PI) *CABS (-K1* (CNl-CN2+i*CN3) -K2* (CN4-CN5+i*CN6) 
+  -T1-T3+CN7+CN8) **2 

SIG=(SIGMA)**.5 
WRITE(16,980)  SIG 
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CONTINUE 

CLOSE(16) 

960 

FORMAT (A) 

970 

FORMAT (15) 

980 

FORMAT (E12. 3) 

STOP 

END 
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APPENDIX  F.  PROGRAM  TO  COMPUTE  a/irb  versus  KB  FOR  SLIT 
CIRCULAR  CYLINDERS 

C  THIS  PROGRAM  COMPUTES  THE  SCATTERING  WIDTH/PI*B  FOR  A  CIRCULAR 
C  CIRCULAR  HOLLOW  CYLINDER  WITH  A  SLIT  OF  ANY  SIZE  FOR  THE 
C  TRANSVERSE  ELECTRIC  POLARIZATION  CASE  USING  THE  ON  -  SURFACE 
C  RADIATION  CONDITION  APPROACH  . 

REAL*8  J(0:200) ,Y(0:200) ,DY(0:200) ,DJ(0:200) ,K1D,K2D 
REAL*8  Jl(0:200) ,Y1(0:200) ,DJ1(0:200) ,DY1(C:200) ,H1 
REAL  EN,NMX,R,U,UU,H,P,PI,PHI,F(99) ,P1 
INTEGER  M,N,NT,NMAX 

REAL  TA,TB,K,K1,K2,TP1,TP2,FIL0N,FR,B,A,L 

COMPLEX  P2,P3,i,Sl,Tl,T2,T3,T4,CBl,CB2,CAl,CA2 

COMPLEX  Cl , C2 , C3 , C4 , C5 , C6 , C7 , C8 , CN2 , CN3 , CN5 , CN6 , CN7 , CN8 

COMPLEX  DN(99) ,S2,S3,S4,S5,CN1,CN4 

EXTERNAL  FI LON 

OPEN (3 , FILE= ' DAT. 71 • ) 

WRITE (*,*)  'ENTER  THE  PHI  ANGLE (REAL)  IN  DEGREES' 

READ(*,*)  P 

WRITE (*,*)  'ENTER  THE  OUTER  RADIUS (REAL)  IN  METERS' 
READ(*,*)  B 

WRITE (*,*)  'ENTER  THE  INNER  RADIUS (REAL)  IN  METERS’ 
READ(*,*)  A 

WRITE (*,*)  '  ENTER  THE  FIRST  SLIT  ANGLE (REAL)  IN  DEGREES' 

READ(*,*)  TA 

WRITE (*,*)  'ENTER  SECOND  SLIT  ANGLE (REAL)  IN  DEGREES' 
READ(*,*)  TB 

WRITE ( * , * )  ' ENTER  NMAX ( INTEGER) ' 

READ(*,*)NMAX 
TPl=TA/57. 29578 
TP2=TB/57. 29578 
PI=3. 1415927 
DTR=PI/180. 

PHI=P*DTR 

i=(0,l) 

M=98 

H=ABS (TP2-TP1) /FLOAT (M) 

C  PLOTTING  ROUTINE 

WRITE (3, 960)  'TE  CASE  PEC  SLIT  CYLINDER  SCATTERING  WIDTH  ' 

NPHI=59 

PHI 1=0.1 

PHI2=6. 

WRITE (3, 970)  NPHI 
WRITE (3, 980)  PHIl 
WRITE (3, 980)  PHI2 
DO  50  R  =0.1,6. , .1 
K1=R*B 
K2=R*A 
K1D=DBLE(K1) 

K2D=DBLE(K2) 

CB1= (8 . *Kl**2-8 . *i*Kl) / (8 . *K1**2-12 . *i*Kl-3 . ) 

CA1=- (8 . *K2**2-8 . *i*K2)/ (8. *K2**2-12 . *i*K2-3 . ) 


56 


CB2=-4 ./ (8 . *K1**2-12 . *i*Kl-3 . ) 

CA2=-4 ./ (8 . *K2**2-12 . *i*K2-3 . ) 

CALL  BES(NMAX,K1D,J,Y,DJ,DY) 

CALL  BES(NMAX,K2D,J1,Y1,DJ1,DY1) 

IF(P.EQ.O.)  THEN 
Tl=i*K* (A-B) *COS (TP2) 

T3=i*K* (A-B) *COS (TPl) 

ELSE 

T1=C0S (TP2 ) / (COS (TP2 ) -COS (TP2+PHI) ) * (EXP ( i*K2* (COS (TP2 ) -COS (TP2+ 
+  PHI) ) )-EXP(i*Kl*(COS(TP2)-COS(TP2+PHI) ) ) ) 


T3=COS (TPl) / (COS (TPl) -COS (TPl+PHI) ) * (EXP(i*K2* (COS (TPl) -COS (TP1+ 
+  PHI) )) -EXP (i*Kl*( COS (TPl) -COS (TPl+PHI) ) ) ) 

ENDIF 

T2=EXP( (-i) *K2*COS (TP2+PHI) ) -EXP( (-i) *Kl*COS (TP2+PHI) ) 

T4=EXP( (-i) *K2*COS (TPl+PHI) ) -EXP( (-i) *Kl*COS (TPl+PHI) ) 

C  EVALUATING  1ST  FUNCTION  F(U) 

DO  100  U  =1,M+1 

UU=TP1+(H*(U-1.) ) 

F(U)=COS(Kl*COS(UU+PHI) ) 

100  CONTINUE 

SI =(0,0) 

P1=0. 

DO  200  N=1,NMAX+1 

NMX=FL0AT(N)-1. 

IF(NMX.EQ.0.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

PI  =  FILON(F,NMX, TPl, TP2,M+1, 99,1) 

WRITE (*,*)  NMX,P1 
P2=EN*i**NMX*J (NMX) 

H1=J(NMX) 

WRITE (*,*)  P2 
P3=P1*P2 
WRITE (*,*)  P3 
SI  =  SI  +  P3 
200  CONTINUE 

C1=S1 

C  EVALUATING  2ND  FUNCTION  F(U)  (CN2=Cl-jC2)  &  CN4 
DO  250  U=1,M+1 

UU=TP1+(H*(U-1.) ) 

F(U)=SIN(Kl*COS(UU+PKI) ) 
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250  CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO  300  N=1,NMAX+1 

NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

P2=(0,0) 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P1=FIL0N(F,NMX,TP1,TP2,M+1,99, 1) 

S3= (-1 . ) **NMX*J1 (NMX) *DJ1 (NMX) 

S4=(-l. ) **NMX*i**(-l. ) *CA1*DJ1 (NMX) **2/ (l.+CA2*NMX**2) 
DN(N) =2 . *PI*EN* (S3-S4) *COS (NMX*PHI) 

P2=EN*i**NMX*J (NMX) 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

300  CONTINUE 

CN4=S2 

C2=S1 

CN2=Cl-i*C2 

C  EVALUATING  3RD  FUNCTION  F(U)  &  CN7 
DO  350  U=1,M+1 

UU=TP1+(H*(U-1.) ) 

F (U) =COS (UU+PHI) *COS (Kl*COS (UU+PHI) ) 

350  CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO  400  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

S5=(0,0) 

P2=(0,0) 

IF(NMX.EQ.0.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

Pl=FILON (F, NMX, TP1,TP2,M+1, 99,1) 

P2=EN*CBl*i** (-NMX-1 . )*(-!.) **NMX*DJ (NMX) / ( 1 . +CB2  *NMX**2 ) 
S3=COS (NMX*TP2 ) *T2 

S4=(CBl*i**(-NMX-l. )*DJ(NMX) )/(l.+CB2*NMX**2) 
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S5= (CAl*l** ( -NMX-1 . ) *DJ1 (NMX) ) / ( 1 . +CA2 *NMX**2 ) 

DN (N) =EN/2 . * (-1. ) **NMX* (S4+S5) *S3 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

400  CONTINUE 

C3=S1 

CN7=S2 

C  EVALUATING  4RTH  FUNCTION  (CN3*C3-jC4)  6  CN8 
DO  450  U=1,M+1 

UU=TP1+(H*(U-1.) ) 

F (U) =COS (UU+PHI) *SIN (Kl*COS (UU+PHI) ) 

450  CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO  500  N=1,NMAX+1 
NMX=FLOAT (N) -1 . 

S3=(0,0) 

S4=(0,0) 

S5=(0,0) 

P2=(0,0) 

IF(NMX.EQ.0.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P1=FIL0N(F,NMX,TP1,TP2,M+1,99,1) 

P2=EN*CBl*i** (-NMX-1 . )*(-!.) **NMX*DJ (NMX) / ( 1 . +CB2*NMX**2 ) 
P3=P1*P2 

S3=C0S (NMX*TP1) *T4 

S4=(CBl*i**(-NMX-l.)*DJ(NMX) )/ (l.+CB2*NMX**2) 

S5= (CAl*i** ( -NMX-1 . ) *DJ1 (NMX) ) / ( 1 . +CA2*NMX**2 ) 

DN (N) =EN/2 . * ( -1 . ) **NMX* (S4+S5) *S3 

S1=S1+P3 

S2=S2+DN(N) 

500  CONTINUE 

CN8=S2 
C4=S1 

CN3=C3-i*C4 

C  EVALUATING  5TH  FUNCTION 

DO  550  U-1,M+1 

UU=TP1+(H*(U-1.) ) 

F(U)=COS(K2*COS (UU+PHI) ) 

550  CONTINUE 

S1=(0,0) 

P1=0. 

DO  600  N=1,NMAX+1 
NMX=FLOAT(N)-l. 
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IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

P1=FIL0N ( F , NMX , TPl , TP2 , M+1 , 99 , 1) 
P2=EN*i**NMX*Jl (NMX) 

P3=P1*P2 
S1=S1+P3 
600  CONTINUE 

C5=S1 

C  EVALUATING  6TH  FUNCTION  (CN5=C5-jC6) 

DO  650  U=1,M+1 

UU=TP1+(H* (U-1.)) 
F(U)=SIN(K2*COS(UU+PHI) ) 

650  CONTINUE 

S1=(0,0) 

P1=0. 

DO  700  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

IF(NMX.EQ.O. )  THEN 
EN=1. 

ELSE 

EN=2. 

ENDIF 

P2=(0,0) 

Pl=FILON (F, NMX, TP1,TP2,M+1, 99,1) 
P2=EN*i**NMX*Jl (NMX) 

P3=P1*P2 
S1=S1+P3 
700  CONTINUE 

C6  ^S1 

CN5=C5-i*C6 


C  EVALUATING  THE  7TH  FUNCTION  &  CNl 
DO  750  U=1,M+1 

UU=TP1+(H*(U-1.)) 

F (U) =COS (UU+PHI) *COS (K2*COS (UU+PHI) ) 
750  CONTINUE 

S1=(0,0) 

S2=(0,0) 

P1=0. 

DO  800  N=1,NMAX+1 
NMX=FLOAT(N)-l. 

S3=(0,0) 

S4=(0,0) 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 

EN=2. 
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ENDIF 

P2=(0,0) 

P1=FIL0N(F,NMX,TP1,TP2,M+1,99,1) 

P2=EN*CAl*i** (-NMX-1. ) * (-1. ) **NMX*DJ1 (NMX) / ( 1 . +CA2*NMX**2 ) 
S3=J(NMX)*DJ(NMX) 

S4=i** (-1) *CB1*DJ (NMX) **2/ (1+CB2*NMX**2) 

DN (N) =2 . *PI*EN* (S3+S4) *COS (NMX*PHI) 

P3=P1*P2 

S1=S1+P3 

S2=S2+DN(N) 

800  CONTINUE 

C7=S1 
CN1=S2 

EVALUATING  THE  8TH  FUNCTION  (CN6=C7-jC8) 

DO  850  U=1,M+1 

UU=TP1+(H*(U-1.) ) 

F(U) =COS (UU+PHI) *SIN(K2*COS (UU+PHI) ) 

CONTINUE 
S1=(0,0) 

P1=0. 

DO  900  N=1,NMAX+1 
NMX=FIiOAT  (N) -1 . 

IF(NMX.EQ.O.)  THEN 
EN=1. 

ELSE 
EN=2. 

ENDIF 
P2=(0,0) 

Pl=FILON (F, NMX, TP1,TP2,M+1, 99,1) 

P2=EN*CAl*i** (-NMX-1 .)*(“!•) **NMX*DJ1 (NMX) / ( 1 . +CA2 *NMX* * 2 ) 
P3=P1*P2 
S1=S1+P3 
CONTINUE 
C8=S1 

CN6=C7-i*C8 

SIGMA=1 ./ (4 . *K1*PI) *CABS (-K1* (CNl-CN2+i*CN3) -K2* (CN4-CN5+i*CN6) 

+  -T1-T3+CN7+CN8) **2 

WRITE (3, 980)  SIGMA 
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CONTINUE 

CLOSE (3) 

960 

FORMAT (A) 

970 

FORMAT (15) 

980 

FORMAT (E12. 3) 

STOP 

END 

C 


850 


900 
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APPENDIX  6.  PROGRAM  TO  COMPUTE  THE  BESSEL  FUNCTIONS 


SUBROUTINE  BES (N, X, J , Y , DJ , DY) 

THIS  SUBROUTINE  WAS  OBTAINED  DURING  EC  4600  FROM  PROFESSOR 
MICHAEL  A.  MORGAN. 

Double  precision  calculation  of  ordinary  Bessel  functions, 
Jn(X) 

and  Yn(X) ,  and  their  first  derivative,  DJ  and  DY,  for  integer 

order  "n”  from  n=0  to  N  with  real  argument  X. 

REAL*8  J(0;200) ,Y(0:200) ,DJ(0:200) ,DY(0;200) , SCALE , JTEMP2 , X 
REAL*8  SCLFAC ,A,B,C,D,E,F,PI, JTEMP , JTEMPl 
PI=3 . 14159265359D0 
IF  (X.EQ.O.ODOO)  THEN 
X  =  0.0  BOUNDARY  CASE 
IF  (N.EQ.l)  THEN 

J(l)  =  O.ODOO 
DJ(1)  =  0.5D00 

ELSE 

DO  5,  I  =  N,  2,  -1 
J(I)  =  O.ODOO 
DJ(I)  =  O.ODOO 
5  CONTINUE 

J(l)  =  O.ODOO 
DJ(1)  =  0.5D00 
ENDIF 

J(0)  =  l.ODOO 
DJ(0)  =  O.ODOO 
Y(N)  =  -l.OD-300 
DY(N)  =  1.0D300 
ELSEIF  (N.EQ.O)  THEN 

POLYNOMIAL  EXPANSION  ONLY  FOR  N  =  0 
CALL  BES0(X, J, Y,PI,DJ,DY) 

ELSE 

RECURSION  FOR  ALL  OTHER  CASES 
Y  IS  A  FORWARD  RECURSION 
CALL  BES0(X, J, Y,PI,DJ,DY) 

Y(l)  =  -DY(0) 

DY(1)=Y(0)  -  Y(l)/X 
IF  (N.EQ.l)  GO  TO  20 
DO  10,  I  =  0,  N-2 

Y(I+2)  =  (2.0D00*(I+1)/X)*Y(I+1)  -  Y(I) 

DY(I+2)  =  Y(I+1)  -  ( (I+2)/X)*Y(I+2) 

10  CONTINUE 

J  IS  A  REVERSE  RECURSION  BASED  ON  A  PAIR  OF  BESSEL  FUNCTION 
POINTS  DERIVED  FROM  A  TRUNCATED  POWER  SERIES  EXPANSION.  THE 
RECURSION  IS  THEN  SCALED  TO  A  KNOWN  VALUE,  J1(X). 

20  SCALE  =  -DJ(0) 

NSAVE  =  N 
IF  (X.LE.N)  THEN 
N  =  5*N+50 
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GOTO  25 

ENDIF 

N  =  IDNINT(N  +  X*X  +  0.5D00) 

C 

25  A  =  1.0DOO/DFLOAT(N+1) 

B  =  1.0D00/DFLOAT(N+2) 

C  =  1.0D00/DFLOAT(N+3) 

D  =  1.0D00/DFLOAT(N+4) 

E  =  1.0D00/DFLOAT(N+5) 

F  =  X/2.0D00 
C 

JTEMP  =  1-A*F**2+0.5D00*A*B*F**4-(1.0D00/6.0D00) *A*B*C*F**6+ 
( 1 .  ODOO/24  .  ODOO)  *A*B*C*D*F**8-  ( 1 .  ODOO/120 .  ODOO)  *A*^-*C*D*E*F**10 
N  =  N  -  1 

A  =  1.0DOO/DFLOAT(N+1) 

B  =  1.0D00/DFLOAT(N+2) 

C  =  1.0D00/DFLOAT(N+3) 

D  =  1.0D00/DFLOAT(N+4) 

E  =  1.0D00/DFLOAT(N+5) 

F  =  X/2.0D00 

JTEMPl  =  1-A*F**2+0.5D00*A*B*F**4-(1.0D00/6.0D00)*A*B*C*F**6+ 
( 1 . ODOO/24 . ODOO) *A*B*C*D*F**8- ( 1 . ODOO/120 . ODOO) *A*B*C*D*E*F**10 
DO  30,  I  =  N+1,2,-1 

JTEMP2  =  2* ((I  -  1)/X)*JTEMP1  -  JTEMP 
IF(DABS(JTEMP2) .GE.1.0D250)  THEN 
JTEMP2  =  JTEMP2*1.0D-250 
JTEMPl  =  JTEMPl*!. OD-250 

ENDIF 

JTEMP  =  JTEMPl 
JTEMPl  =  JTEMP2 
IF  ( (1-2) .LE.NSAVE)  THEN 
J(I-2)  =  JTEMP2 

ENDIF 
30  CONTINUE 
C  SCALING 

N  =  NSAVE 

SCLFAC  =  SCALE/J(1) 

DO  40,  I  =  0,  N 
J(I)  =  SCLFAC*J(I) 

IF  (I.EQ.O)  THEN 
GOTO  40 

ENDIF 

IF  (ABS(J(I)/J(I-1) ) .LT.l.OD-50)  THEN 
J(I)  *  J(I)*1.0D250 

ELSEIF  (ABS(J(I)/J(I-1) ) .GT.1.0D50)  THEN 
J(I)  =  J(I)*1.0D-250 

ENDIF 
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DJ(I)  =  J(I-l)  -  (I/X)*J(I) 

40  CONTINUE 
ENDIF 
RETURN 
END 
C 

SUBROUTINE  BES0{X, J,Y,PI,DJ,DY) 

C  FOR  ZERO  ORDER  BESSEL  FUNCTIONS  ONLY 

DIMENSION  J (0:200) ,  Y(0:200),  DJ(0:200),  DY(0:200) 

DOUBLE  PRECISION  J,  Y,  X,  PI,  DJ,  DY,  FO,  FI,  THETAO 
DOUBLE  PRECISION  THETAl,  A 
IF  (X.LE.3.0D00)  THEN 
A  =  X/3.0D00 

J(0)  =  l.ODOO  -  2.2499997D00*(A**2)  + 

1.2656208D00*(A**4)  - 

+0 . 3163866D00* (A**6)  +  0 . 0444479D00* (A**8 ) 

0.0039444D00*(A**10)  + 

+0.00021D00*(A**12) 

Y(0)  =  (2.0D00/PI)*DLOG(X/2.0D00)*J(0)  +  0 . 36746691D00 

+ 

+. 60559366D00* (A**2)  -  0 . 74 3 50 3 8 4D00* ( A**4 )  + 

0.25300117D00*(A**6) 

+  -  0.04261214D00*(A**8)  +  0 . 00427916D00* (A**10) 

0.00024846D00* 

+  (A**12) 

DJ (0) =-X* ( . 5D00-0 . 56249985D00* (A**2) +0 . 21093573D00 
+*(A**4)-  0.03954289D00*(A**6)  +  0. 00443319D00* (A**8) 

0.00031761 

+D00*{A**10)  +  0.00001109D00*(A**12) ) 

D  Y  (  0  ) 

(“1 . ODOO/X) * ( (2 . ODOO/PI) *X*DLOG(X/2DOO) * (-1. ODOO* 

+DJ (0) ) -0 . 6366198DOO+0 . 2212091D00* (A**2) +2 . 1682709D00* (A**4) 

+1 . 3164827D00* (A**6)  +  0 . 3 1 2 3 9 5 1 DO 0 * ( A* * 8 ) 

0.0400976D00*(A**10) 

+  +  0.0027873D00*(A**12) ) 

ELSE 

A  =  3. ODOO/X 

FO  =  .79788456D00  -  0 . 00000077D00*A 

0.00552740D00*(A**2) 

+-0. 00009512D00* (A**3)  +  0 . 0 0 1 3 7 2 3 7 DO 0 * ( A* * 4 ) 

-0 . 00072805D00* (A**5) 

++0.00014476D00*(A**6) 

THETAO  =  X  -  0.78539816D00  -  0 . 04166397D00*A  - 

0.00003954 

+D00*(A**2)  +  0.00262573D00*(A**3)  -  0 . 00054125D00* (A**4 )  - 
+0.00029333D00*(A**5)  +  0 . 00013558D00* (A**6) 

J(0)  =  FO*DCOS (THETAO) /DSQRT(X) 

Y(0)  =  FO*DSIN (THETAO) /DSQRT(X) 

FI  =  0.79788456D0C  +  0 . 00000156D00*A  +  0. 01659667D00*A*A 
++0.00017105DOO*(A**3)  -  0 . 00249511D00* (A**4 )  +  0 . 00113653D00 
+*(A**5)  -0.00020033D00*(A**6) 
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THETAl  =  X  -  2.35619449D00  +  . 12499612D00*A  +  0.00005650 
+D00*(A**2)  -  0.00637879D00*(A**3)  +  0. 00074348D00* (A**4)  + 

+  0.00079824D00*(A**5)  -  0.O0029166D0OO* (A**6) 

DJ(0)  =  -F1*DC0S(THETA1)/DSQRT(X) 

DY(0)  *  -F1*DSIN (THETAl) /DSQRT(X) 

ENDIF 

RETURN 

END 
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APPENDIX  H.  PROGRAM  TO  COMPUTE  THE  INTEGRATION  ROUTINES 


B 

This  program  computes  INTEGRAL  F(x)  cosTx  dx  if  KEY  =  1  and 
B  A 

Int  F  (X)  sin  Tx  dx  when  KEY=2  using  Filon's  mwthod.  F(x)  is 
A 

assumed  to  be  given  as  a  table  of  functional  values  in  the 
array  at  H  equidistant  points  from  A  to  B,  M  odd.  MAX  is  the 
dimension  of  the  array  declared  in  the  calling  program. 

Ref:  P.  J.  Davis  and  P.  Rabinowitz,  Methods  of  Numerical 
Integration,  Academic  Press,  p.  488 


REAL  FUNCTION  FILON  (F,  T,  A,  B,  M,  MAX,  KEY) 

INTEGER  KEY,  M,  N,  I,  MAX 

REAL  F  (MAX),  AL,  BE,  GA,  SUM,  TH,  H,  S,  C,  SI,  Cl,  S2,  C2 

REAL  FI,  F2,  SU,  SUl,  Al,  A,  B,  T,  TH2,  TH4 

N  =  M  -  1 

H  =ABS(B  -  A)  /  FLOAT  (N) 

TH  =  T  *  H 
S  =  SIN  (TH) 

C  =  COS  (TH) 

TH2  =  TH  *  TH 

IF  (ABS  (TH)  .LE.  0.1)  THEN 
TH4  =  TH2  *  TH2 

AL  =  2.  /  45.  *  TH  *  TH2  *  (1.  -  TH2  /  7.  +  TH4  /  105.) 

BE  =2.  /  3.  +2.  /  15.  *  TH2  -  4.  /  105.  *  TH4 

GA  =  4.  /  3.  -  2.  /  15.  *  TH2  +  1.  /  210.  *  TH4 

ELSE 

AL=1.  /TH+S*C/  TH2  -  2.  *  S  *  S  /  (TH2  *  TH) 

BE=(2.  +2.  *C*C  -4.  *S*C/  TH)  /  TH2 

GA  =  4.  *  (S  /  TH  -  C)  /  TH2 

END  IF 

SUM  =0.0 

FI  =  F  (1) 

F2  =  F  (M) 

51  =  SIN  (A  *  T) 

52  =  SIN  (B  *  T) 

Cl  =  COS  (A  *  T) 

C2  =  COS  (B  *  T) 

Al  =  A  +  H 

GO  TO  ( 1 ,  2 ) ,  KEY 
1  SU  =  F2  *  S2  -  FI  *  SI 

SUl  =  -0.5  *  (F2  *  C2  -  FI  *  Cl) 

DO  3  I  =  2,  N,  2 

SUM  =  SUM  +  F  (I)  *  COS  (Al  *  T) 

Al  =  Al  +  H 

SUl  =  SUl  +  F  (I+l)  *  COS  (Al  *  T) 

3  Al  =  Al  +  H 

GO  TO  4 
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If) 


2  IF  (T  .EQ.  0.)  THEN 

FILON  =  0. 

RETURN 
END  IF 

SU  *  FI  *  Cl  -  F2  *  C2 

SUl  =  -0.5  *  (F2  *  S2  -  FI  *  SI) 

DO  5  I  =  2,  N,  2 

SUM  =  SUM  +  F  (I)  *  SIN  (A1  *  T) 

A1  =  A1  +  H 

SUl  =  SUl  +  F  (I+l)  *  SIN  (A1  *  T) 

A1  =  A1  +  H 

FILON  =  H  *  (AL  *  SU  +  BE  *  SUl  +  GA  *  SUM) 
RETURN 
END 
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